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Abstract. We present the results of a numerical investigation of the turbulent kinematic dynamo problem in a high 
Prandtl number regime. The scales of the magnetic turbulence we consider are far smaller than the Kolmogorov 
dissipative scale, so that the magnetic wavepackets evolve in a nearly smooth velocity field. Firstly, we consider 
the Kraichnan-Kazantsev model (KKM) in which the strain matrix is taken to be independent of coordinate 
and Gaussian white in time. To simulate the KKM we use a stochastic Euler-Maruyama method. We test the 
theoretical predictions for t he gr o wth of rates of the m agnetic energy and higher order mo ments llKazantsevi 
1 19681: iKulsrud fc Anderson! Il992t IChertkov et all Il999l). the shape of the energy spectrum jKazantsev ~ ^68t 
lKTnsru^^^AndersorT^)9"a Sc^tochihi^^^il^ ]2002al : iNazarenko et all l2003t) and the behaviour of the polari- 
sation and spectral flatness ( N^iz£a"antc^^lTT200^T !hi general, the results appear to be in good agreement with 
the theory, with the exception that the predicted decay of the polarisation in time is not reproduced well in the 
stochastic numerics. Secondly, in order to study the sensitivity of the KKM predictions to the choice of strain 
statistics, we perform additional simulations for the case of a Gaussian strain with a finite correlation time and also 
for a strain taken from a DNS data set. These experiments are based on non-stochastic schemes, using a timestep 
that is much smaller than the correlation time of the strain. We find that the KKM is generally insensitive to 
the choice of strain statistics and most KKM results, including the decay of the polarisation, are reproduced well. 
The only exception appears to be the flatness whose spectrum is not reproduced in accordance with the KKM 
predictions in these simulations. 
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1. Introduction 

Many astrophysical magneti c fields are believed to be gen- 
erated via a dynamo actio n jMoffattl fl978l IParkeri fl979t 
IChildress fc Cxilberttfl99i in systems where the magnetic 
field diffusivity k is much less than that of the kinematic 
viscosity v. These systems correspond to large magnetic 
Prandtl number Pr = v / k, flows. Such astrophysical situ- 
ations are observed in the interstellar medium as well as in 
protogalactic plasma clouds where the magnetic Prandtl 
num ber varies respectively from Pr ~ 10 14 to Pr ~ 
10 22 llChandrarAll998tlKulsrudLll999l:ISchekochihin et all 
2002a). Such high values of Pr lead to an interesting in- 
terval of scales (from 7 to 11 decades in our examples) 
below the viscous cut-off, but above the magnetic diffu- 
sive scale, where magnetic fluctuations are stretched by a 
randomly changing smooth velocity field. The small-scale 
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kinematic dynamo problem can thus be formulated as to 
whether a small initially "seeded" magnetic field, subject 
to stretching by a prescribed smooth velocity field, will 
grow in time. 

One can picture the evolution of an ensemble of mag- 
netic wavepackets, each traveling together with the fluid 
particles which are distorted by the local strain. By as- 
suming a given form of the strain statistics, an idealised 
homogeneous and isotropic strain that is a Gaussian white 
noise process, one can simplify the problem, leading to a 
productive framework from which analytical results can 
be obtained. This formali sm is known as the Kraichnan- 
Kazantsev model (KKM) l|Kraichnan fc Nagaraianill967t 
iKazantsevt ^968). It should be remembered that such a 
model is a simplification and a real turbulent velocity 
field will exhibit finite correlations of the fluid velocity 
field, as well as being to some degree non-Gaussian, inho- 
mogeneous and anisotropic. Therefore, in this numerical 
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investigation we will also consider more realistic examples 
of the velocity field. 

The work presented here is based on a kinematic 
model, making use of the linear induction equation of the 
magnetic field B(x, t) 



D t B = B • Vv + reV 2 B, 



(1) 



courier space one-point correlators wmcn nave rccenti; 
been investigated analy tically l|Schekochihin et al 1 12002s 
Nazarenko et al , 2003). In the past much emphasis ha 



where D t = dt + v • V is the material derivative, v the 
velocity, and re the magnetic diffusivity. Physically, we 
can have in mind the following picture for the evolution 
of our dynamo system. Initially, for small re, the B-ficld 
will be " frozen" into th e flow behaving as a passive vec- 
tor field ijSowardL Il994j) . After some period of time, and 
hence growth of the magnetic field, diffusion will become 
important and will reduce this rate of growth. In partic- 
ular, using the KKM the growth of th e total mean mag- 
netic energy < B 2 (i) > was obtained bv lKazantsevI lfl968) 
and that of th e highe r order moments < B 2ra (t) > by 
IChertkov et"aH l|l999(L ISchekochihin et alJ (|2002al) have 
extended this work to include the effects of compressibil- 
ity. 

The work presented here has two primary aims. The 
first aim is to n umerically verify the analytical scal- 
ing la ws found by iKazantsevI l)l968|) and IChertkov et alJ 
l[l999l) for the even moments < B 2n (i) > (where n is 
a positive integer) in both the perfect conductor (non- 
diffusive) and diffusive regimes. We will also compute the 
Fourier space one-point corre lators which have recently 

M 

has 

been placed solely on investigations of the magnetic energy 
spectrum. However, higher order correlators are also of 
significance. Indeed, they describe statistics of the fluctua- 
tions around the mean energy distributions in the Fourier 
space. Also, it has been shown that the new quantities, 
the mean polarisation of the magnetic turbulence and the 
spectral flatness, are of particular interest as they char- 
acterise the properties of the small-scale intermittency 
l|Nazarenko et al.ll2003j) . These newer theoretical objects 
will also be investigated numerically. Our second aim is to 
study the sensitivity of the dynamo model to the choice of 
strain statistics. In a quest to test the universality of the 
analytical results found using the KKM, we will consider 
the case of a finite-correlated Gaussian strain field. We will 
also consider a more realistic strain field generated from a 
Direct-Numerical-Simulation (DNS) of the Navier-Stokes 
equations. 

The presentation of this paper has been organised into 
six main sections. In the next section we give a brief in- 
troduction to the small-scale dynamo problem, with a 
summary of the KK M a nalytical results a s derived by 
IChertkov et~aT] l|l999l) and lKazantsevI (|l968h . In the third 
section we outline some of the new KKM analytical re- 
sults in regard to the behaviour of the one-point Fourier 
space correl ators; a more deta i led di scussion of which can 
be found in INazarenko et all l|2003|) . The fourth section 
provides an overview of the numerical method used to 
investigate the KKM and its subsequent results. In the 



fifth section we outline the method and results of the 
finite-correlated Gaussian and DNS based investigations. 
Finally, we draw our conclusions in the sixth section. 

2. Moments of B 

Let us briefly summarise the analytical results for 
the moments o f B(£ ) obtained within the KKM by 
IChertkov et al ] Jl999l). At scales below the viscous cut- 
off Batchelor l|Batcheloil 11950. 1954) argued the velocity 
field would be random and smooth (Batchelor regime). In 
MHD a uniform velocity cannot alter a system's magnetic 
field. The flow 



(2) 



where a(t) is a coordinate independent random strain ma- 
trix, represents one of the simplest choices of a velocity 
field that allo ws for the transfer of kin etic energy into mag- 
netic energy l|Zel'dovich et al.Lll98# . The incompressibil- 
ity condition V • v = in this case corresponds to Tr[a] = 
for all time. It is natural to reformulate our mathemati- 
cal description in terms of Lagrangian dynamics (namely 
following particles ) in which case D+ = d / dt in the induc- 
tion equation (TjQ) . IZel'dovich et al.l l)l984^ llSowa,rdLlT994^ 
showed that given an initial condition, the induction equa- 
tion Q can be re-written in Fourier space as 



B(M) = J(i)B(q,0) exp 
where the Jacobian J(t) satisfies 
a J, 



k 2 (t')dt' 



dt 



(3) 



(4) 



with the initial condition J(0) = /, with I the unit matrix, 
and q = k(0) is the initial wavevector related to k(t) by 

q = J T k. (5) 

Moments of B are calculated via two independent averag- 
ings, one over initial statistics, an d the other over variou s 
realisations of the strain matrix (|Chertkov et all Il999h . 
The initial B field is assumed to be homogeneous, isotropic 
and Gaussian, with zero mean and the following variance 



<B a (k,Q)B b (k',0) > = 
k a k b \ , 

Oab 7TT- I k tin 



(6) 



S(k - k') 



where L is the length-scale at which the "seed" magnetic 
noise is initially concentrated and E m is a constant which 
determines the initial tur bulence intensity . After averaging 
over the initial statistics IChertkov et al.l l)l999|) obtained 
the following expression for B 2 (t) 

B 2 (i) = j E m e-o 2L2 e~ 2 ^ A V J ia ($ ab - S^j J tt dq, (7) 

with summation over repeated indices and where q = |q| 
and A(t) satisfies 

A(i) = / dt' J- l (t')J- T {t'). (8) 
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The full even moments < B 2 ™(i) > where calculated by 
taking the value of B 2 to the correct power and then av- 
eraging over the realisations of the velocity held. For an 
isotropic, delta-correlated, Gaussian strain field 



< <J ab (t)a ab (0) > = lOAiJ(t), 



(9) 



where Ai is the growth rate of a material line element, cor- 
responding to t h e larg est Lyapunov exponent of the flow, 
IChertkov et alJ l)l999(l derived the following expressions 
for the scaling behaviour of the even moments < B 2 "(i) >. 
In the perfect conductor regime t < 3td/ (2n + 3), we have 



< B 2 ™(<) > 



2Ain(2n+3)t/3 



(10) 



where on dimensional grounds the scale at which mag- 
netic diffusion becomes important is r& ~ yj k/X% and the 
corresponding time is tj ~ X[ \xi{L/r^). In the diffusive 
regime t > 3td/(n + 2), we have 



< B 2 "(i) > 



Am(n+4)t/4 



(11) 



For n = 1, this result was derived earlier bv iKazantsevI 
(1961). 



3. Fourier Space Correlators 

Traditionally, Fourier space analysis of the dynamo 
problem has only really involved investigations of 
the fc-space correlator correspondin g to the e nergy 

spectrum E(k,t) = < IB Otll 2 > (IKazantsevI 11968 

Kulsrud fc Andersonl Il992t ISchekochihin et all l2002a|) 
Kaz antsevl (|1968) studied the dynamo system as an 



eigenvalue problem from which he was able to obtain 
the growth exponents of the total magnetic energy. 
The evolution of the energy spectrum in fc-space has 
also been studied both an alytically and numerically by 
iKulsrnd fc Andersonl JTflflj. 

In this KKM analysis, the way in which we define the 
incompressible strain matrix is different to used in 
IChertkov et"afl |l99<J. We write 



^1% 

3 3 



where G is a 3 x 3 matrix made up of independent 
Gaussian N(0; 1) random variables, normally distributed 
(i.e. Gaussian) with zero mean and a standard deviation 
of 1, such that 



< G l3 (t)G a {0) >=S tk 5 jl S{t). 



(13) 



This choice is motivated by its convenience in numerical 
modelling. One should note there is a degree of flexibility 
in how we define our strain statistics, a more detailed dis- 
cussion of which can be found in the appendix. The KKM 
results remain the same for different choices of the strain 
statistics, the only difference being that time is re-scaled 
by a constant factor. 



3.1. Energy Spectrum 
3.1.1. Energy Spectrum 



Solution 



In the limit of zero diff usion k = 0, the 3D e nergy spec- 
trum has the solution (Nazarenko et all 120031) 



-1/2 



exp 



_31n !M )A /50A 

Ant J \ 4 / 



,(14) 



where q is the wavenumber at which E(k) is initially con- 
centrated at t = 0, and Eq is a constant depending on 
the initial conditions. If we let td be the diffu sion time 
as defined in section 121 or IChertkov et af] l|l999f) . then for 
large time Qr 1 lrv 1 {k/q) -C t -C id, i- c - m the perfect con- 
ductor regime, we have E ~ fc -1 / 2 . Alternatively, we can 
think of the fc" 1 / 2 slope as being present over an interval 
of wavenumber space fc c _ <C k <ti k c+ where 



k c ±(t) ~ exp(±-v/4fti/3). 



(15) 



The critical wavenumbers k c+ and /c c _ govern when the 
exponential log squared term above becomes important. 
We see the fronts at k c+ and k c ^ propagate exponentially 
in time, and within this fc -1 / 2 range the spectrum also 
grows exponentially. By integrating equation (|14f) . over 
the whole of wavenumber space, one can obtain the growth 
rate of the total magnetic energy 



< 



B 2 (t) >= J E(k,t)dk~ const exp(^j^). 



(16) 



which agrees with result (|10fl for n = 1 when we take into 
account the slightly different definitions of the strain co- 
variance. Indeed, one finds that the X±t of (|10|) and fit of 
(|16|l are equivalent if w e re-scale time by a fac tor of 4/5. 
Hence the scaling laws of lChertkov et all ljl99fl|) , presented 
in sectional for the diffusive regime should be re- written 
for our choice of strain statistics as 



< B 2n (t) > ~ g2n O (2n+3)t/3 



Further, in the diffusive regime we have 



(12) < B 2n (t) > 



Qn(n+4)t/4 



(17) 



(18) 



3.1.2. Energy Spectrum : k/0 Solution 

If all wavevectors are initially of the same length 
' q|, and we are in t erested in large time a symp totics 



1 Schekochihin et all l2002al iNazarenko et all l2003h . the 

solution in the diffusive regime takes the form 



E(k,t) 



const 



ffito/4 k -V 2 t- 3 '*Ko(k/k K ), 



(19) 



where Kq is a MacDonald function of zeroth order and k K 
is a constant given by k K = yCtf&K. In fact this k K is the 
wavenumber at which diffusion becomes important. It is 
therefore equivalent to the estimated kd ~ l/ra derived 
on purely dimensional grounds in section [3 
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3.2. Flatness and Polarisation 

Although the energy spectrum solutions described above 
provides us with some very useful information, they do not 
provide us wit h a complete picture of o ur system. In the 
recent work of iNazarenko et all l|2003f) higher one-point 
correlators of the magnetic field in Fourier space were 
studied 1 . In this paper, as well as for the energy spectrum 
E(k,i), we will also study two other correlators 



S(k,t) =< |B(k)| 4 >, 
T(k,t) =< |B 2 (k)| 2 > 



(20) 
(21) 



Let u s briefly review the results from INazarenko et alJ 
l)2003h regarding these correlators. 

3.2.1. Polarisation Spectrum and Flatness : k = 
Solutions 

In the perfect conductor regime, the fourth order 
correlators of Fourier amplitudes S =< |B(k)| 4 > and 
T =< |B 2 (k)| 2 > have exact solutions 



S(k,t) 



T(k,t) 



i n 
Vt U 

Vbexp 
1 (k 

vi U 

Vo exp 



1/2 

exp 



3 In 



\k/q) \ 



4m j 



(22) 



7 



Q exp 



mt \ 



1/2 



exp 



31n 2 (fc/g) \ 



(23) 



4 

yQoexp 



mt\ 



where Vo and Qq are constants. In large time the Qq terms 
can be neglected and these spectra develop a k 1 / 2 scaling 
between two propagating fronts fc c _ <C k <C k c+ . 

Combining these two fourth order correlators, we find 
an important new quantity P = (S — T) / S , the physical 
significance of which becomes clear when we write P as 



P 



< IBI 4 > 



(24) 



< IBI 4 > 



\B 4 



where <pi is the phase of the complex component Bi and 
the operator takes the imaginary part of an expres- 
sion. In this form we see that P contains information about 
the phases of the fc-space modes. P can be thought of as 
the mean normalised polarisation. Indeed, P = corre- 
sponds to the plane polarisation of the Fourier modes. In 



1 These were objects of the form < |B| 2 "|B 2 | 2m > which 
represent a basis for all one-point correlators in magnetic tur- 
bulence which is isotropic. 



contrast for a Gaussian magnetic field one finds the polar- 
isation P = 1/3. Combining (|22|l and (|23|l . in large time 
we find the normalised polarisation P evolves as 



P= ^exp(-6fii). 
Vo 



(25) 



That is, the polarisation is independent of wavenumber k 
and decays exponentially. This solutions tells us that in 
the perfect conductor regime, the Fourier modes of the 
magnetic field will eventually become plane polarised. In 
comparison, we should recall a Gaussian field has a finite 
polarisation, thus our turbulent magnetic field is far from 
being Gaussian. 

Another important measure of turbulent intermittency 
in a fluid, both in real and Fourier space, is the flatness. 
The spectral flatness is defined as the ratio 



F 



F? 



< IB(k) 



> 



< |B(k)P >2 



(26) 



Using H26|) . one finds a Gaussian field has a constant flat- 
ness of 3/2. On the other hand small-scale intermittency 
corresponds to a field with a flatness that grows both in 
time and in wavenumber space. Indeed, in the perfect con- 
ductor regime, the mag netic field displays such s mall-scale 
intermittent behaviour l(Nazarenko et al.ll2003h . with 



F~ t 



1/2 



exp 



p/ 2 



exp 



3 In (k/q) \ 

Ant J 



(27) 



For fc c _ < k < k c+: we see there is a range of fc 3 / 2 scaling, 
while for larger k there is an increase in the flatness arising 
from the log-squared exponential term. This intermittency 
in Fourier space can be attributed to the presence of co- 
herent structures, as will be discussed later. 



3.2.2. Polarisation Spectrum and Flatness : n / 
Solutions 

In the diffusive regime the fourth order correlators S and 
T have large time solutions l|Nazarenko et all l2003j) 



S ~ k 1 ' 2 ^' 2 
T ~ k 1 ' 2 ^' 2 



y 0e 21O t /4 + £ Qoe -30 t /4 



Vo e 2im ^ - ^Q e- 3Ot/4 



k Rl : 



.(29) 



where fc K .2 = \JQ/12k. This fc K .2 can again be interpreted 
as the wavenumber at which diffusion becomes important. 
However, in comparison to the energy spectrum we should 
note that the spectral cut-off will be at a smaller wavenum- 
ber as k K .2 < k K . Importantly we see that the normalised 
polarisation (|25|l behaves identically in both the diffusive 
and perfect conductor regimes. In contrast, the behaviour 
of the flatness is modified 



F = 



_S_ 

E 2 



t,3/2.3/2-iim Ko(k/k Ki2 ) 
(K (k/k K )Y 



(30) 
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when compared to l|27|) . For small k, we again have a re- 
gion of k z / 2 scaling but now with a log correction aris- 
ing from the MacDonald functions. For large k, below the 
spectral cut-off, the additional MacDonald functions act 
to heighten the flatness. That is, the introduction of a 
finite diffusivity actually increases the small-scale inter- 
mittency. 



4. KKM Numerical Investigation 

The key to modelling the KKM numerically is the suc- 
cessful integration of equation Q to find J. As the strain 
matrix on the right hand side of contains noise, it 
is evident that this is not a normal ordinary differential 
equation (ODE) and should instead be interpreted in the 
stochastic sense. Following the analytical formalism (|12|l 
and (|13|l . the elements of the strain matrix a will be made 
up from a matrix of Gaussian random variables G, such 
that the incompressibility condition Tr[er] = is satis- 
fied. As with any stochastic formalism one must decide 
the form of the integral. Mathematically both the Ito and 
Stratonovich definitions of the stochastic integral are cor- 
rect. It is only in a few special cases that both interpreta- 
tions produce the same solutions. In the present problem 
this is not the case. Indeed, the Stratonovich formalism is 
the correct interpretation for our problem as we are us- 
ing white noise as an i dealisation of what i n reality is a 
coloured noise process l(Kloeden et all Il997^ . 

The simplest numerical scheme one can use to solve 
Stochastic Differential Equations (SDE) is called the 
Euler-Maruyama scheme (EMS), and it can be thought 
of as generalisation o f the normal forward Euler method 
l|Kloeden et all ^997). The stochastic numerical experi- 
ments presented here can be divided into two separate 
implementations. Both these codes require the stochas- 
tic integration of the J equation Q), which we perform 
numerically using an EMS. Averaged quantities are cal- 
culated over many realisations of the strain statistics and 
initial conditions. An ensemble of particles is chosen as the 
initial condition, randomly distributed on a unit sphere in 
wavenumber space |q| = 1 (corresponding to an isotropic 
distribution in real space). Each particle is subjected to 
its own realisation of the strain matrix, and its wavevector 
evolves according to equation J^J. 

4.1. Code 1 

Using the integration process described above for the evo- 
lution matrix J, Code 1 (CI) determines how the magnetic 
field evolves according to equation . As each particle is 
subjected to its own realisation of the strain matrix the 
integral Q over initial conditions simplifies to just the 
value of the integrand. Higher order quantities, such as 
B 4 ■ ■ ■ , can then be calculated before performing the final 
averaging over strain statistics. The initial magnetic noise 
is determined by © with L = 1 . In the case of k ^ 
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one must also determine the value of the matrix A(t) from 
equation JSJ. This is achieved via a forward Euler scheme. 

A jk (t n+1 ) = A jk (t n ) + At Jil l {tn)Ju T {tn)^ ■ (31) 

We see therefore, that we must also determine the value 
of J -1 . This can be done in two different ways. Firstly, 
we could calculate the inverse of J at each step in the 
usual manner by finding the adjoint and determinant of 
J. However, as some elements of J can become very large 
in time one needs to be careful in finding the determinant 
as round-off errors can be easily introduced. Alternatively, 
we can integrate the separate stochastic equation for the 
evolution of J -1 

^{J- 1 ) = -(J->, (32) 

using an EMS. 

The primary aim of this code is to verify the scaling 
laws for the even moments < B 2n (t) >, (|17|) and (|18|) . 
However, by recording the values of the integrand of (7J 
and the corresponding wavenumbers k = |k| for each par- 
ticle, we can also construct snapshots of the energy spec- 
trum in time. 

CI has also been set up to calculate the eigenvalues 
of the matrix J J T so that the Lyapunov exponents can 
be calculated. Indeed, one important result for systems 
with random time-dependent strain is that, for nearly 
all realisations of the strain, the matrix (lit) ln( J T J) is 
found to stabili s e at l arge time IColdhirsc h et all Il987t 
iFalkovich et all 1200 lj) . The eigenvectors U in this case 
tend to a fixed orthonormal basis for each realisation, 
where i = 1,2,3. As the Jacobian J at t = is the iden- 
tity matrix, J T J will also take the form of the identity 
initially. In time the strain a will distort J via (0J) . An ini- 
tial sphere will be deformed into an ellipsoid, the volume of 
which will be conserved by the incompressibility condition 
Det[J] = 1. The length of the principal axis of the ellipsoid 
correspond to the eigenvalues e 2pi of J T J, while the eigen- 
vectors fi give the directions of these axis. It is evident that 
in considering a time-dependent stochastic strain we have 
moved to a system corresponding to Lagrangian chaos. In 
which case, we can define the Ly apunov exponents Xj in 
terms of the l i miting eigenvalues jGoldhirsch et all . Il987t 
ISowardl Il994t IFalkovich et all l200l|) as 

A,= lim lln(ff(J T J)fi), (33) 

t^oc it 

for i = 1,2, 3. It is customary to order the Lyapunov expo- 
nents in terms of size, Xi > A2 > A3. A material element 
aligned along one of the eigenvectors will in the asymp- 
totic limit expand or contract at the rate exp(Aii). The 
Lyapunov exponents do not depend on the realisations of 
the strain, while in contrast the asymptotic eigenvectors 
are realisation dependent. Incompressibility ensures that 
Y^ = i = which in turn tells us that one Lyapunov 
exponent must be positive if two or more are non-zero. A 
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positive Lyapunov exponent corresponds to exponential 
growth of a material line element. 

In the non-degenerate case, if the statistics of the 
strain matrix are symmetric to time reversal, which is 
true if the strain is assumed delta c orrelated as is the 
case for the KKM then A 2 = l)Kraichnanl 11974 
Balk ovskv fe Fouxonl. fl999) and he nce the incompressibil- 
ity condition gives A3 = — Ai < 0. IChertkov et alJ l)l999() 
found that it is the largest Lyapunov exponent Ai that is 
responsible for the g r owth of the moments < B 2n (i) >. 
Bal kovskv fe Fouxonl l|l999|) have also considered the ma- 
trix JJ T , which evolves in time according to 



dt 



(JJ T ) = a(JJ T ) + (JJ T )a 7 



(34) 



This differs to the matrix J T J in that it does not stabilise 
to some fixed axis at large time. Indeed, the eigenvectors 
will continue to rotate for every realisation of the strain 
matrix. As in . the Lyapunov exponents can also be 
found from the logarithm of the eigenvalues of JJ T at 
large time. The eigenvalues of JJ T and J T J, and hence 
their Lyapunov exponents, are the same as they share the 
same characteristic polynomial. 

To determine the Lyapunov exponents of our system 
numerically we need to determine either JJ T or J T J at 
each timestep. This can be achieved by either calculating 
JJ T via its own dynamic equation (|34|l and integrate it 
using an EMS, or alternatively calculate JJ T or J T J at 
each timestep from J. In practice either approach works 
equally well. Next one must determine the eigenvalues of 
the symmetric matrix J J T which must be real. Although 
the matrix is only 3x3 solving the characteristic polyno- 
mial using the solution to a cubic equation is impractical 
as the growth/reduction of elements J in time produce 
round-off errors numerically. This in turn leads to com- 
plex eigenvalues and further inaccuracies. 

The method employed in this numerical study is that 
of Jacobi transformations IjPress et al l ll993h . The key to 
this method is a sequence of similarity transformations 
(rotations). In the context of ellipsoids, this algorithm can 
be thought of as a method by which the axis of the sys- 
tem are rotated and re-aligned to lie along the principle 
axis of the ellipsoid, the length of these axis correspond- 
ing to the eigenvalues, and the directions to the eigen- 
vectors. Although the Jacobi rotation method is not per- 
fect, it does however perform better than the previously 
mentioned cubic equation approach, in that the eigenval- 
ues remain real for all time. Given the set of eigenvalues 
exp(2pi) of JJ T , we define the time-dependent Lyapunov 
exponents to be 



Xi(t) = — < ln(exp(2 j o 4 )) >= - < Pl >, 



(35) 



where i — 1,2, 3. Analytically the Lyapunov exponents 
are independent of the strain realisations at large time, 
but the averaging above has been performed over all the 
realisations and initial conditions to improve accuracy. 



4.1.1. k = simulations 

We will start by setting n = in the numerics so that we 
can investigate the perfect conductor regime. Firstly, let 
us compare the scaling of < B 2 (i) > with theoretical so- 
lution (|17J) . The left-hand graph of Fig. ^ corresponds to 
the n = 1 case and we see the numerically generated slope 
agrees nicely with the theoretical scaling (dashed line). 
For the higher order moments < B 2n (i) >, good agree- 
ment with the analytical scaling (| 1 T|l can also be seen. 
The right-hand graph in Fig.^shows the slope for n = 2. 
As all moments are calculated from the integral J7J, it 
should be noted that any fluctuations in B 2 (t) will get am- 
plified at higher orders, consequently the resulting graphs 
will become progressively noisier. The only solution to this 
problem is to increase the number of the realisations over 
which averaging is performed. 

Let us now investigate the energy distribution in 
wavenumber space. The energy spectrum < B 2 (fc,i) > is 
obtained at a particular instant in time by recording each 
particles wavenumber k = |k| and "weight" corresponding 
to the integrand of 10. One then constructs a histogram 
by first finding the largest and smallest wavenumbers k max 
and k m i n of the all the particles, and then dividing the in- 
terval log k m in < log k < log k m ax into a finite number 
of bins. Particles are then sorted into these bins and their 
corresponding weights are summed. Finally, we normalise 
the total weight in each bin, by dividing through by its bin 
width. One must also divide through by a factor of 47rfc 2 
which originates from the solid angle integration required 
in conve rting a ID spectr um into that of a 3D one (see for 
example iMcCombl l|l990l) L 

In Fig. |3 the energy spectrum has been plotted for two 
different times t = 9 (left figure) and t — 18 (right figure). 
Here 57 = 0.16 and thus using the definition ljl5|l for the 
critical wavenumber fc c +(i) we obtain log 10 k c ±(9) — ±0.6 
and log 10 fc c ±(18) = ±0.85. In Fig. the same spectrum 
has been plotted for t = 27 (left figure) and t = 31.5 
(right figure), correspondingly log 10 k c+ (27) = 1.04 and 
lo 

S10 &c+(31.5) — 1.13. One should recall that theoret- 
ically for wavenumbers fc c _ <C k <C fc c + we expect a 
kr x l 2 scaling range 2 , and this is indeed the case here 
(the straightline in each graph corresponding to a slope 
of fc -1 / 2 ). As a whole we can see that the spectrum moves 
vertically upwards in time. This corresponds to the time 
dependent increase of the magnetic energy at fixed fc 
predicted in the theoretical result <|14(1 . This large time 
asymptotic solution has also been plotted on each figure 
and we see good agreement with the numerically gener- 
ated histograms. 

It should be noted that the later figures show the spec- 
trum for log 10 fc > 0. This is because the front propa- 
gating to small fc becomes extremely noisy at large time. 
The reason for this will become apparent shortly. It is 
evident from these results that it would be beneficial if 



Here, we ar e con sidering a 3D spectrum. 
iKulsrud fc Anderson! il992T) considered the equivalent 
ID spectrum which has a corresponding fc 3 ' 2 scaling regime. 
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we could extend the number of decades over which this 
scaling is observed. A possible approach would be to run 
the simulations for a much longer time, remembering that 
logfc c ±(t) ~ \/fti. In doing this one has to be careful of 
numerical stability and the generation of round-off errors 
as some elements of J will get ever larger. However, this is 
not the only problem, we should remember that our dis- 
tribution of particles, initially concentrated at our initial 
wavenumber |q| = 1, will broaden in time. The peak of 
this distribution will move towards larger k, each individ- 
ual particle performing its own random walk in wavenum- 
ber space. Fig. 0] shows the distribution of particles at 
four successive times t = 9, 18, 27 and 35 for a simula- 
tion with Q = 0.16. The graph has been normalised so 
that the area under the curve is one. We see that with 
log wavenumber coordinate, the particle distribution fits 
a moving Gaussian profile. The mean of this profile moves 
to higher wavenumbers as fit, while the standard devi- 
ation grows as Regrettably, the energy spectrum's 
k~ x / 2 scaling range of interest lies in one of the tails of 
the particle distribution. This means that there will be 
fewer particles available for averaging when calculating 
the energy spectrum histogram in this region. 

The only way to improve this averaging is to increase 
the total number of particles used in a simulation. Of 
course, this is not an efficient way of improving the statis- 
tics, as increasing the number of rarer particles located 
in each tail, will correspond to a much greater increase in 
the number of particles located in the centre of the particle 
distribution. This is an inherent problem with this type of 
particle simulation. Also, as the particle distribution gets 
more spread out in time, fluctuations at both ends of the 
energy spectrum histogram will increase. This effect can 
clearly be seen in Fig. and That is why the front fc c - 
that propagates to smaller k < q < 1 in time becomes so 
noisy, as the bulk of the particles are travelling to larger 
wavenumbers. For this reason in the remaining results we 
will concentrate only on the log 10 k > region of fc-space. 



4.1.2. K ^ simulations 

We will now introduce a non-zero n into the numerical 
model to investigate the effects of diffusion. The first thing 
we will investigate here is the scaling behaviour of the total 
magnetic energy < B 2 (t) > in the diffusive regime. The 
left-hand graph of Fig.[S]shows the results of a simulation 
with £1 = 0.36 and k = 0.005. The two separate scaling 
regimes are clearly apparent and the analytical slopes (|17f) 
and l|18fl . for the perfect conductor and diffusive regimes 
respectively, are in good agreement with the numerics. It 
is inevitable that to generate similarly smooth graphs for 
the higher order moments, one would have to use a greater 
number of realisations. FigureEJand the right-hand graph 
of Fig.[5]show the scaling of the next three even moments. 
As expected these graphs are noisier, but still in agreement 
with the analytical results. 



Next we will consider the energy spectrum. The in- 
troduction of a finite diffusivity will cut-off the energy 
spectrum at some finite wavenumber. Figure [7J shows two 
successive snapshots of the energy spectrum at times t = 6 
and 12 for a simulation with k = 0.005 and = 0.36. As 
expected, one observes a spectral cut-off and for the results 
shown in Fig. d the parameter choice gives k K = yJVt/Qn. 
Theory predicts that for k < k K there should be a re- 
gion of fc -1 / 2 scaling (|19|l (with a logarithmic correction). 
For this particular simulation we have log 10 k K ~ 1.1. As 
with the perfect conductor regime, the spectrum is seen 
to move vertically upwards corresponding to a steady in- 
crease in the total magnetic energy. The large time the- 
oretical solution lj*H?|l has also been plotted in Fig. [7J and 
good agreement is observed. 

We will briefly review the Lyapunov exponents gener- 
ated in these CI numerical experiments because it help 
us to understand the origins of some numerical errors in 
our method. Using the Jacobi transform method described 
earlier, one can calculate the eigenvalues e 2pi from either 
J T J or JJ T , and then order them in terms of size so 
that pi > p2 > P3- The Lyapunov exponents are then 
determined via (|35|l . Figure [S] shows < 2pi{t) > for a sim- 
ulation with fl = 0.36. As expected the graph tends to a 
flat line in time, and Ai > 0. We would expect Ai ~ ft, 
that is < lne 2pi >= 2£Ai = 2tCfl where C is the time- 
scaling constant mentioned earlier. In Fig. [SJa line corre- 
sponding to C — 4/5 has also been plotted. From theory 
we would expect the remaining two Lyapunov exponents 
to take the values A2 = and A3 = — Ai, with conser- 
vation of volume corresponding to Ai = 0. The left 
and right-hand graphs of Fig. El shows < 2p 2 (t) > and 

< 2ps(t) > respectively. Looking at the range of the verti- 
cal axis we see A2 does indeed stay around zero but there 
is a gradual deviation from the theoretical value in time. 
On the other hand, A3 is not the same as — Ai. The line 

< hie 2 '' 3 >= -2tAi = -2tCfl with C = 4/5 has also been 
plotted for comparison. Numerically, we see that ^ Ai 7^ 
and hence, according to the Lyapunov exponents, the vol- 
ume is not conserved. It should also be mentioned that in 
our simulations, the determinant of the Jacobian Det[J], 
is held within 1% of unity, which means that the volume 
is well conserved. We conclude therefore that the devia- 
tions observed in the numerical values of Ai from theory 
are mainly due to errors generated in finding the eigen- 
values of JJ T or J T J using the Jacobi transform method. 
This is not overly surprising if we remember that some 
elements of JJ T are very large in comparison to others. 
Indeed, that is why the alternative method of finding the 
eigenvalues via solving the cubic characteristic polynomial 
of JJ T runs into difficulties. In the context of the ellip- 
soids, essentially the problem arises because each ellipsoid 
is becoming ever larger in one direction and smaller in 
another. 

It should be noted that the magnetic Reynolds num- 
ber is relatively small in these simulations. The magnetic 
Reynolds number is by definition R m = ^/L 2 f2/ n, giv- 
ing a range of values 32 to 72 in the above simulations, 
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where we have taken L to be of unity. Thus, it is only 
the relative sizes of VL and k that are important. Higher 
Reynold number simulations are possible but again in- 
creasing Q will make the stochastic fluctuations greater 
and hence one will need to reduce the timestep in each 
simulation. Alternatively, decreasing the diffusivity k will 
delay the time at which the second regime commences. In 
which case, the corresponding simulation will have to be 
run over a larger time interval. 

4.2. Code 2 

Code 2 (C2) uses the same foundations as CI, an ensemble 
of particles are initially distributed on the surface of a unit 
sphere in wavenumber space. However, now each particle 
is given a randomly orientated magnetic field satisfying 
B • q = 0. That is, the real and imaginary parts of the 
magnetic field are randomly orientated in a plane which 
is at right-angles to the particle's initial wavevector and 
tangent to the unit sphere in fc-space. Particles are again 
allowed to evolve in wavenumber space via but in 
contrast to CI, information about the full magnetic field 
is retained by the integration of equation 0. In the case 
of k ^ we need to determine the integral in © . As with 
CI this is achieved through the use of a forward Euler 
scheme. In this way one can investigate the spectra of 
Fourier space quantities of interest such as the magnetic 
energy, polarisation and flatness. 

4.2.1. k = simulations 

Starting with the energy spectrum < B 2 (fc, t) >. This has 
been plotted in the left-hand graph of Fig. ^] along with 
the large time theoretical curve l|14|l and slope fc -1 / 2 , for 
the choice of parameters = 0.16 and t = 35. As with 
the CI experiment Fig. the numerical and analytical 
results support each other well. As with the investigation 
of higher order moments < B 2 (i) > in section 14.1.11 the 
higher order correlators, such as S =< |B| 4 (fc,t) >, will 
generate noisier histograms in comparison to their lower 
order counterparts. Nevertheless, the following numerical 
results are still very much of interest. Firstly, we will ex- 
amine the correlator S which has been plotted as the right- 
hand graph of Fig. ^] The large time asymptotic solutions 
has also been plotted, and there appears to be good quali- 
tative agreement with the numerics, although some devia- 
tion is evident at large k. As with the energy spectrum the 
analytical results tell us that there exists a fc +1 / 2 scaling 
range for k < k c+ . In this simulation log 10 k c+ = 1.13. A 
slope of fc +1 / 2 has been included in the graph and again 
there would seem to be agreement. However, it must be 
said that the fc-interval for this scaling range is perhaps 
too short and the histogram too noisy to be conclusive. 

For the same simulation, Fig.llllshows the normalised 
mean polarisation P (left figure) and spectral flatness 
F (right figure). Let us consider the P spectrum. From 
the theoretical results we expect P to be independent 



of wavenumber. Although, the histogram is noisy, there 
does indeed appear to be no /c-dependence, and the re- 
sulting spectra is flat. However, for the KKM we also 
have the prediction that the polarisation tends to zero in 
time (I25|l . Unfortunately, this is not the case numerically. 
Indeed, taking an average over all k in the histogram we 
can find an approximate value for P from the above spec- 
trum. Repeating this procedure for various snapshots in 
time we find the normalised polarisation remains steady in 
time at the value of 0.32. One should remember that for a 
Gaussian field analytically we found P = 1/3. Therefore, 
in our numerical simulation the polarisation seems to be 
that of a Gaussian field. However, the corresponding spec- 
tral flatness is far from being Gaussian. Indeed, the flat- 
ness F in this perfect conductor simulation can be seen in 
the right-hand graph of Fig. El Below log 10 k c+ — 1.13 we 
see agreement with the analytical solution of k +3 / 2 . While 
at larger k the theoretical increase in the flatness, due to 
the log squared term in (|27|l . is also well represented. 

One may well ask, why is there a discrepancy between 
theory and the numerically generated polarisation, but not 
the flatness? It helps if we can understand the physical sig- 
nificance of these two quantities. We return to the familiar 
example of the evolution of an initial ball of isotropic mag- 
netic turbulence is wavenumber space. For each realisation 
this ball is deformed into an ellipsoid with one large, one 
short and one neutral direction. One can visualise this el- 
lipsoid as an elongated flat cactus leaf with thorns aligned 
to the magnetic field direction. Note that in this picture 
one component of the magnetic field (transverse to the 
cactus leaf) is dominant which is captured by the fact 
that the polarisation P introduced above, tends to zero at 
large time. Another consequence of this picture is that the 
wavenumber space will be covered by the ellipsoids more 
sparsely at large k. This produces large intermittent fluc- 
tuations of the Fourier transformed magnetic field. These 
fluctuations are quantified by the growth of the flatness F 
as fc 3//2 , a clear indication of this small-scale intermittency. 

Numerically, it would appear that the flatness is a little 
more robust than the polarisation, arising naturally in the 
numerics as the ellipsoids become increasingly sparse at 
large k. In contrast, the apparently Gaussianity of the 
polarisation can be attributed to the mis-alignment of the 
cactus needles. Indeed, the magnetic needles should align 
themselves with the smallest principal axis of the ellipsoid. 
In this direction the cactus leaf is becoming continuously 
thinner in time, and hence numerically one might expect 
some errors to creep into the alignment of B. Indeed, this 
seems to be related to the errors observed in calculating 
the smallest Lyapunov A3. 

In fact, although on average our numerical simulations 
appear to mis-represent the behaviour of the polarisation 
completely, things are not as bad as they would initially 
seem. Indeed, further insight can be obtained by consid- 
ering the behaviour of an ensemble of particles that are 
subjected to the same realisation of the strain. Figure IT2l 
shows the real magnetic fields of a set of 500 wavepackets 
that have been subjected to two different realisations of 
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the strain matrix 3 . In the left-hand figure it is clear the 
magnetic field in this realisation is far from being plane 
polarised, the magnetic vectors still being predominantly 
random in their orientation at this given point in time 
t = 6. In contrast, in the right-hand figure, which has 
also been taken at t = 6, we see the ellipsoid has become 
very elongated and the magnetic field appears plane po- 
larised. Earlier snapshots, at t = 2 and 5, of the magnetic 
field for this particular realisation can be found in Fig. 
1131 and the evolution of the corresponding polarisation 
spectra is shown in Fig. ^] In this case, we see the po- 
larisation does indeed decay in time and is far from being 
Gaussian by t = 6. In some respects, therefore, our numer- 
ical model does get some of the polarisation's behaviour 
at least qualitatively correct. 

Before we consider the effects of finite diffusivity, we 
should check how the growth rates of the spectra behave 
in time. That is, how rapidly they move vertically up or 
down in time. By recording the value of each spectra for 
k = 1, at different snapshots in time, we can compare 
this numerical growth to the theoretical solutions l|14l) . 
123, 123) and 123 for E, S, T and F respectively. These 
results can be found in Fig. Along with the numeri- 
cal values (points) the corresponding theoretical growth 
has also been plotted (solid line). It should be noted that 
this is only a rough check. Firstly, because we have had 
to ignore the \fi prefactor in each analytic solution, and 
secondly because the numerical values of the spectra at 
k — 1 are likely to fluctuate at large time as we are again 
in the tails of the particle distribution. Nevertheless, all 
the spectra are growing in time in agreement with theory. 

4.2.2. k 7^ simulations 

In this section we will discuss the results of a set of C2 
simulations with a finite diffusivity of n = 0.005. FigurelTCI 
shows snapshots of < B 2 (k, i) > at the times t — 6 and 12. 
As with the CI simulations, the finite diffusivity produces 
a spectral cut-off. As expected the spectrum moves verti- 
cally upwards in time, with the second histogram getting 
noiser as the density of particles in the interval k < k c + is 
reduced. Comparing these energy spectra to their equiva- 
lent CI counterparts, Fig. [7| we should note that the C2 
numerics appear to be slightly hyper-diffusive. This effect 
is most probably numerical, originating from the integra- 
tion of the k 2 integral found in © ■ 

The higher order correlators S and T have been plot- 
ted in Fig. 1171 for a time t = 0. We should note that the 
particles here have been sorted into 100 bins. However, 
the axis in each graph has been re-scaled to take into ac- 
count the spectral cut-off. For small k, qualitatively these 
figures are in agreement with the large time asymptotic 
solutions. However, as with the energy spectrum, these 
results seem to slightly hyper-diffusive. The final figure 
we will consider is 1181 which shows flatness F at t — 6. 

3 The imaginary magnetic fields are qualitatively similar and 
have therefore not been included. 
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Again we observe good qualitative agreement with the- 
ory. In particular, one can observe a k 3 ^ 2 scaling region for 
small k and a heightened flatness at larger k in agreement 
with equation l|30jl • This coincides with the analytical pre- 
diction that the inclusion of a finite diffusivity increases 
the flatness (and hence small-scale intermittency) at large 
k. The corresponding normalised mean polarisation spec- 
trum has not been included here as it is again flat, much 
like the perfect conductor regime (see Fig. 111(1 but with 
a spectral cut-off. As with the non-diffusive case, on av- 
erage the normalised polarisation appears Gaussian with 
P = 0.32. Again, it is interesting to consider the magnetic 
field of an ensemble of particles that have been subjected 
to the same realisation of the strain matrix. Figure HTA 
shows two such snapshots of an ensemble of 500 magnetic 
particles. These two figures should be compared with Fig. 
I12l which show the same realisations but without diffusion. 
The right-hand figure in particular demonstrates why, in 
the diffusive regime, an ellipsoid will cover fc-space more 
sparsely due to the decay of its magnetic field at its tips. 
This is the reason why spectral flatness increases in the 
diffusive regime. 

5. Numerical Investigations Beyond the KKM 

The purpose of the stochastic codes outlined above was to 
investigate the small-scale dynamo system for the case of a 
Gaussian white strain field. In contrast, the following two 
codes were developed to test the universality of the ana- 
lytical KKM results summarised in section [21 when more 
realistic representations of the velocity field are chosen. 
These simulations are based on the integration of the in- 
duction equation Q , which when written in a Lagrangian 
frame of reference in fc-space 4 takes the form 

-^B rn = a ml Bi - nk 2 B m , (36) 

where d/dt = dt + kjdkj and 

kj = &ijki. (37) 

Numerically we consider an ensemble of particles 
(wavepackets), whose individual magnetic fields evolve ac- 
cording to i|36|) and wavevectors according to (|37l) . 

5.1. Code 3 

In code 3 (C3), we consider a synthetic Gaussian strain 
field with a finite correlation time. The r.m.s. of the dif- 
ferent strain components in this numerical experiment 
ranged from 2.9 to 3.6 and the correlation time was 

4 Strictly speaking, the Fourier transforms used in deriving 
this equation must be taken over a box. The box size is greater 
than the characteristic length-scale of magnetic turbulence, but 
less than the characteristic length-scale of the velocity field. In 
this case, model |(2J corresponds to ignoring the quadratic (in 
x) terms in the scale separation parameter. The strain here is 
measured along the fluid path. 
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0.02. Thus the correlation time was about sixteen times 
less than the characteristic strain distortion time. If we 
choose a timestep which is much smaller than the corre- 
lation time, we no longer need to consider a stochastic 
scheme, but can instead use a more traditional second- 
order Runge-Kutta method to integrate and i|37|) . 
This experiment consisted of 2048 strain realisations with 
2048 magnetic wavepackets and the timestep was cho- 
sen to be twenty times smaller than the strain correla- 
tion time. Initially, the particles are randomly distributed 
within a ball of radius 2 in the /c-space. Each particle 
is given a randomly orientated magnetic field (with ran- 
dom phase) that is transverse to its initial wavevectors q 
and which has a random amplitude less than or equal to 
2 x 10~ 4 . 

Figure [20| shows the energy spectrum (left figure) and 
the growth of the total energy (right figure). The en- 
ergy spectrum is consistent with theory, although there 
is no clear A; -1 / 2 range. At later time (in the diffusive 
regime) one can see that the spectrum acquires a cut- 
off. In the right-hand figure, the total energy grows expo- 
nentially with two different scaling regimes, as expected. 
The time at which diffusion becomes important is around 
4. In agreement with theory, in the diffusive regime the 
growth rate of the total energy is approximately equal to 
the growth of the individual fc-modes. Remarkably, the ra- 
tio of the growth rate in the perfect conductor regime to 
the growth rate in the diffusive regime is in good agree- 
ment with theoretical predictions. In Fig. [2] one can see 
the corresponding mean polarisation (left figure) and spec- 
tral flatness (right-figure) from this simulation. One can 
see that the polarisation spectrum is flat and it rapidly 
decreases in time. Both of these features are in agree- 
ment with the theoretical KKM results. Thus, the KKM 
based analytical prediction that the magnetic turbulence 
becomes plane polarised is much better represented in this 
simulation than in the previously discussed stochastic sim- 
ulations. The spectral flatness rapidly grows in time which 
is also in agreement with the KKM theory. However, there 
is no obvious region of the theoretical fc 3 / 2 scaling. 

5.2. Code 4 

Code 4 (C4) is almost identical in implementation to C3. 
However, in this numerical experiment the strain matrix 
components were obtained from a 512 3 spectral DNS of 
the Navier-Stokes equations at a Taylor Reynolds number 
of approximately 200. The strain time series calculated via 
recording 



along 512 fluid paths. At each of these fluid particles we 
put 4096 magnetic wavepackets. The r.m.s. of the differ- 
ent strain components in this experiment ranged from 6.6 
to 9.3 and the correlation time was approximately 0.08. 
Thus, the correlation time in this case is of the same or- 
der as the inverse strain rate which is natural for real 



Navier-Stokes turbulence where both values are of the or- 
der of the eddy turnover time at the Kolmogorov scale. 
The timestep in this simulation was taken to be two hun- 
dred times smaller than the strain correlation time. 

In Fig. 123 we consider the energy spectrum (left fig- 
ure) and growth of the total energy (right figure) . There 
is a good agreement with the theoretical behaviour of the 
energy spectrum. One can see a region of k 1 / 2 scaling and 
a log correction in the diffusive regime. Note that the pure 
k 1 / 2 scaling changes to the log-corrected spectrum at the 
time when the high-fc tail hits the dissipative scale which 
approximately corresponds to the third curve in Fig. 1221 
measured at t = 0.44. The growth of the total energy is 
again well represented. The growth rate exponent changes 
at around t = 0.4 which marks the transition between the 
perfect conductor and the diffusive regimes. In agreement 
with the theory, the growth of the total energy the diffu- 
sive regime is approximately equal to the growth of the 
energy found in individual k modes. Again the ratio of 
the growth rates in the perfect conductor and the diffu- 
sive regimes is consistent with theoretical predictions. 

The C4 results for the mean polarisation and the flat- 
ness spectra are shown in Fig. [23 One can see that the 
mean polarisation (left figure) decays in time and is spec- 
trally flat which is in agreement with the KKM theory. 
The spectral flatness (right figure) grows in time, but the 
shape of the spectrum is different from the theoretical re- 
sult of fc 3 / 2 . This is similar to the corresponding flatness 
results obtained in C3 simulations. In fact, the deviations 
of the flatness from theory, in both the C3 and C4 simu- 
lations, might not be due to changes in strain statistics. 
The correlation time of the strain in the C3 simulations 
is short, and therefore one would expect that its numer- 
ical results would be closer in behaviour to KKM than 
C4. However, the degree to which the flatness differs from 
the KKM results is similar in both the C3 and C4 simula- 
tions. Thus, we may conclude that our algorithm is failing 
to reproduce some features of the higher order correlators. 

6. Discussion 

In this numerical investigation we have used a range of 
numerical models, based on Lagrangian particle dynam- 
ics, to investigate the small-scale turbulent dynamo prob- 
lem. For the KKM based stochastic models, using a delta- 
correlated Gaussian strain, we firstly considered the scal- 
ing of the even order moments < H 2n (t) > and found 
good agreement with the theoretical predictions (|17|) and 
(|18fl for the perfect conductor and diffusive regimes re- 
spectively. Next we considered various spectral objects in- 
cluding the magnetic energy, polarisation and flatness, and 
made comparisons to their large time asymptotic analyt- 
ical solutions. Because of the rich analytical results avail- 
able, the KKM is a good testing ground for our stochastic 
numerical models. The behaviour of the energy spectrum 
in both diffusive and non-diffusive simulations was found 
to grow in time and fit the large-time asymptotic forms 
predicted by theory, with evidence of a k^ 1 / 2 scaling re- 
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gion in each case. The fourth order correlators were also 
considered. Although noisier than the energy spectrum, 
these spectra also appeared consistent with the analyt- 
ical results. However, there was evidence of some hyper- 
diffusion in these numerical results are high wavenumbers. 
The new analytical objects, the mean polarisation and 
spectral flatness, were investigated. The flatness was well 
behaved with an observed k 3 ^ 2 scaling and the growth 
rate in time in agreement with the analytical predictions. 
In contrast, the polarisation of the magnetic turbulence 
in these simulations was mis-represented. Although it was 
spectrally flat, in agreement with the theory, it appeared 
to saturate at its Gaussian value of 1/3 rather than de- 
caying in time. However, on considering individual real- 
isations it was apparent that at least some of the pre- 
dicted decaying behaviour was observed qualitatively in 
these simulations. This also provided us with a useful vi- 
sualisation of the dynamo system in terms of a cactus leaf 
in fc-space. The linear polarisation corresponds to the k- 
space alignment of the magnetic field "thorns" along the 
direction of the smallest Lyapunov exponent, i.e. trans- 
verse to the cactus leaf pad. The flatness growth with k 
is directly linked to the fact that at larger k the cactus 
leaves will cover the fc-space more sparsely. We also com- 
puted the Lyapunov exponents and found that the largest 
exponent was well represented which ensures the correct 
growth of the total magnetic energy. The behaviour of the 
smallest exponent was badly reproduced, which may be 
connected to the observed behaviour of the polarisation. 

A second set of experiments was also performed to 
test the universality of the KKM based theoretical re- 
sults with respect to changes in the strain statistics. These 
two simulations, used a finite-correlated Gaussian and a 
DNS generated strain fields respectively. For the finite- 
correlated Gaussian strain based model, the energy spec- 
trum was found to be consistent with theory, although no 
fc~ 1//2 scaling region could clearly be distinguished. For 
the DNS strain based model there was an agreement with 
theory for both the growth rate of the total energy and 
the fc -1 / 2 scaling. Remarkably, these two numerical exper- 
iments demonstrated excellent agreement with the theo- 
retical behaviour of the polarisation, in each case being 
both spectrally flat and decaying in time. Also in agree- 
ment with the theory is the rapid time growth of the flat- 
ness. In contrast, the spectral shape of the flatness in these 
simulations was not as well represented as in the stochas- 
tic based results. For both the finite-correlated Gaussian 
and DNS strain based simulations the flatness spectrum 
was far from the predicted k 3 ! 2 scaling. We believe this 
difference from theory is a numerical artifact from our al- 
gorithm, rather than directly due to changes in the form 
of the velocity field. We conclude that the gross features 
predicted by the KKM theory seem to be robust and in- 
sensitive to chan ges in the strain stat is tics. 

Interestingly. ISchekochihin et all l)200ll l2002b|) re- 
cently investigated the behaviour of the real space cur- 
vature of the field C = b • Vb where b = B/|B|. They 
found that the curvature of the magnetic field is anti- 
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correlated with its strength, which corresponds to folded 
and strongly stretched structures. This agrees with our 
result that the Fourier modes of the magnetic field tend 
to a state of plane polarisation. However, it should be 
noted that the Fourier polarisation gives more information 
than the curvature statistics. Indeed, the zero curvature 
requirement allows both layered and filamentary configu- 
rations. Further, layered magnetic field may vary its direc- 
tion as one passes from one layer to another. Our analysis 
narrows down the choice and indicates that the magnetic 
field structures are layers in coordinate space, thus ruling 
out the existence of any filament structures. Further, the 
plane polarization result also eliminates the possibility of 
any "twists" in the magnetic field between the individual 
layers. That is, the sheet magnetic field is aligned along 
the same direction. Although the field in successive sheets 
can be parallel or anti-parallel to each other. In fact, the 
presence of one neutral direction in the Lagrangian defor- 
mations tells us that the layers have a finite width in one 
direction and thus look, in coordinate space, like ribbons 
with the magnetic field lying along these ribbons. 

Acknowledgements. We would like to thank Warwick 
University's Fluid Dynamic Research Centre for the use of 
their computer facilities. 

Apppendix 

In general, the strain matrix can be defined as 

< o-«(t)(r H (0) >= (39) 

Ci(S ik 5ji + C 2 SijS k i - (1 + C 2 )Sii5j k )5(t), 

where Ci are constants (see for example iMcCombl ()l990|) 'l. 
Written in this way the strain statistics automatically sat- 
isfy the incompressibility condition. 

For our particular choice (|12|l . we have, 

< o-«(t)<r w (0) >= \SikSjt - *(*)> (40) 

giving C\ — fl and Ci = —1/3. In cont rast, the covari- 
ance c hosen bv lChertkov et all l)l999(l and lFalkovich et all 

E3 is 

< <Tij(t)a kl (Q) >= Y^ik^jl ~ Sijhl ~ SuS jk )S(t), (41) 

corresponding to C\ = 4Ai/3 and Ci = —1/4, where Ai 
is the growth rate of a material line element in the flow. 
Setting i = k = a and j = I = /3 in both (|40Jl and l|41(l we 
find, 

< o-«,j(tK/s(0) >= 8flS(t), (42) 

< <r a p{t)cr a p{Q) >= 10Ai«5(t). (43) 

The latter ex pressi on is the previous stated 
IChertkov et alJ l)l999h formalism from section 
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Fig. 1. log < B 2 (t) > (left figure) and log < B 4 (i) > 
(right figure) for a CI simulation, with SI = 0.16 and k = 
0. Averaging was performed over 600000 realisations. The 
theoretical solutions l(T7|) have also been plotted (dashed 
lines) . 
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Fig. 2. E =< B 2 (fc,t) > for a CI simulation with SI = 
0.16 and k — 0. Averaging has been performed over 600000 
realisations. The time is t = 9 (left figure) and t = 18 
(right figure). The particles have been sorted into 100 bins. 
The theoretical solution (|14(l (curve) and a slope of A; -1 ' 2 
(straight line) have also been plotted. 
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Fig. 3. E =< B 2 (fc,i) > for a CI simulation with SI = 
0.16 and k — 0. Averaging has been performed over 600000 
realisations. The time is t = 27 (left figure) and t = 31.5 
(right figure). The particles have been sorted into 100 bins. 
The theoretical solution Ijl4(l (curve) and a slope of k~ x /' 2 
(straight line) have also been plotted. 
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Fig. 4. Normalised particle distribution N(k,t) for a CI 
simulation with f2 = 0.16. Averaging has been performed 
over 120000 realisations. Snapshots are shown at t = 9 
(top left), t = 18 (top right), t = 27 (bottom left) and 
t = 35 (bottom right). The particles have been sorted into 
100 bins. A Gaussian curve has also been fitted to the 
distribution. 





Fig. 5. log < B 2 (t) > (left figure) and log < B 4 (i) > 
(right figure) for a CI simulation, with £1 = 0.36, k = 
0.005. Averaging was performed over 120000 realisations 
(left figure) and 480000 realisations (right figure). The 
theoretical slopes (|17|l and (|18fl have also been plotted 
(dashed lines). 




Fig. 6. log < B 6 (t) > (left figure) and log < B 8 (t) > 
(right figure) for a CI simulation with SI = 0.36 and 
k = 0.005. Averaging has been performed over 480000 re- 
alisations. The theoretical slopes l|17[) and l|18fl have also 
been plotted (straight lines). 
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Fig. 7. E =< B 2 (fc,i) > for a CI simulation with ft = 
0.36 and k — 0.005. Averaging has been performed over 
480000 realisations. The time here is t — 6 (left figure) 
and t = 12 (right figure). The particles have been sorted 
into 100 bins. The theoretical solution l|19|) (curve) and a 
slope of fc -1 / 2 (straight line) have also been plotted. 




Fig. 8. < 2pi(t) > for a CI simulation with fl = 0.36. 
Averaging has been performed over 480000 realisations. A 
slope of (8/5)flt has also been plotted (dashed line). 
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Fig. 9. < 2p 2 {t) > (left figure) and < 2p 3 (t) > (right 
figure) for a CI simulation with £1 = 0.36. Averaging 
has been performed over 480000 realisations. A slope of 
— (8/5)0* (right arrow) has also been plotted (dashed 
line) . 
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Fig. 11. P=(S-T)/S (left figure) and F = S/E 2 (right 
figure) for a C2 simulation with £1 = 0.46 and k = 0. The 
time here is t — 35. Averaging has been performed over 
600000 realisations. The particles have been sorted into 
400 bins. In the right figure the theoretical curve 1)2 7[) and 
slope fc +3 / 2 (straight line) has also been plotted. 
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Fig. 10. E =< B 2 (k,t) > (left figure) and S =< 
|B| 4 (fc, t) > (left figure) for a C2 simulation with Q = 0.46 
and k — 0. The time here is t — 35. Averaging has been 
performed over 600000 realisations. The particles have 
been sorted into 100 bins. The theoretical solution (|14|l 
(curve) and a slope of fc -1 ' 2 (straight line) have also been 
plotted in the left figure. While the curve (|22|) and a slope 
of fc +1 / 2 (straight line) have also been plotted in the right 
figure. 




Fig. 12. The magnetic field of 500 wavepackets in a C2 
simulation with il = 0.36 and k — 0. The figures show 
the particle's positions in /c-space at t = 6 and their cor- 
responding real magnetic fields. The left and right figures 
correspond to different realisations of the strain matrix. 
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Fig. 13. The magnetic field of 500 wavepackets in a C2 
simulation with Q — 0.36 and k = 0. The figures show 
the particle's positions in /c-space and their corresponding 
real magnetic fields at t = 2 (left figure) and t = 5 (right 
figure) for one realisation of the strain matrix. The strain 
field is the same here as the realisation used in right-hand 
graph figure El 
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Fig. 14. The polarisation P = Q/S of 500 wavepack- 
ets in a C2 simulation with = 0.36 and k = 0. Here 
all the particles are subjected to the same realisation of 
the strain matrix. The graph shows the polarisation at 
t = 2 (top line), t = 5 (middle line) and t — 6 (bottom 
line). The corresponding snapshots of the particle's posi- 
tions and magnetic field can be found in figures H3l and 
respectively. 
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Fig. 15. Growth of E(k, t) (top left), F(k, t) (top right), 
5(k, t) (bottom left) and T(k, t) (bottom right) in time 
at k = 1 for a C2 simulation with = 0.36. The theo- 
retical growth for E, F, S and T, (23, El and (HJ 
respectively, have also been plotted (solid lines). 
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Fig. 16. E =< B 2 (fc, t) > for a C2 simulation with Q = 
0.36 and k = 0.005. Averaging has been performed over 
480000 realisations. The time here is t — 6 (left figure) 
and t = 12 (right figure). The particles have been sorted 
into 100 bins. The theoretical solution i|19|) (curve) and a 
slope of k~ x / 2 (straight line) have also been plotted. 
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Fig. 17. S =< |B| 4 (M) > (left figure) and T =< 
|B 2 | 2 (fc,i) > (right figure) for a C2 simulation with Q = 
0.36 and k = 0.005. The time here is t = 6. Averaging 
has been performed over 480000 realisations. The parti- 
cles have been sorted into 100 bins, (but the axes have 
been re-scaled). The theoretical curves (I28|l (left figure) 
and l(2*S|) (right figure) as well as a slope of fc +1 / 2 (straight 
line) have also been plotted. 




Fig. 19. The magnetic field of 500 wavepackets in a C2 
simulation with CI = 0.36 and k = 0.005. The figures show 
the particle's positions in fc-space at t = 6 and their corre- 
sponding real magnetic fields for two different realisations 
of the strain matrix. The left figure corresponds to the 
same strain field as the left-hand graph of figure ^| and 
similarly for the right-hand graphs. 
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Fig. 18. F = S/E 2 for a C2 simulation with Q = 0.36, 
K = 0.005 and t = 6. Averaging has been performed over 
480000 realisations. The particles have been sorted into 
100 bins, (but the axis have been re-scaled). The theoret- 
ical curve (|30() and slope k +3 ^ 2 (straight line) have also 
been plotted. 
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Fig. 20. The magnetic energy spectrum E(k,t) (left fig- 
ure) and growth of the total magnetic energy < B 2 (t) > 
(right figure) for the case of a synthetic Gaussian strain 
field with finite correlation time. The right-hand figure 
also shows the growth of the energy for a selection of in- 
dividual fc-modes. 
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Fig. 21. The normalised mean polarisation (left figure) 
and spectral flatness (right figure) for the case of a syn- 
thetic Gaussian strain field with finite correlation time. 



Energy 




Fig. 22. The magnetic energy spectrum E(k, t) (left fig- 
ure) and growth of the total magnetic energy < B 2 (i) > 
(right figure) for the case of a strain obtained from a 512 3 
DNS. The right-hand figure also shows the growth of the 
energy for a selection of individual fc-modes. 




Fig. 23. The normalised mean polarisation (left figure) 
and spectral flatness (right figure) for the case of a strain 
obtained from a 512 3 DNS. 



